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Abstract 

The distortion of a liquid droplet from its spherical shape has a profound influence on its drag coefficient. The areas 
of greatest uncertainty are the drag at high Reynolds and Weber numbers as well as acceleration and history effects. 
Though several reduced models have been developed for use in spray modeling, a detailed, direct numerical simula- 
tion offers the opportunity to learn more about the underlying physics. A numerical investigation has been per- 
formed to study the variation of the droplet drag with change in deformation factor of the liquid drop, and to better 
understand the dependence on Ohnesorge number and acceleration. The transient nature of these simulations will 
advance the understanding of flow effects on the droplet drag. 


Introduction 

The accuracy of any spray model depends greatly 
on the correctness of the drag coefficient computed for 
the liquid droplets. Experiments have shown that drop- 
lets having a Weber number significantly larger than 
unity tend to deform, which results in a change in their 
drag characteristics. The distortion of the droplet causes 
the increase in the frontal projected area, causing a de- 
viation from the ideal spherical drag characteristics. 

Until recently droplet distortion and its effect on 
drag characteristics have been studied experimentally 
using wind tunnels [1] and shock tube experiments [2]. 
A summary of these experiments and the correlations 
obtained by them can be found in [3]. These experi- 
ments are limited by the range of droplet sizes, fluid 
densities and viscosities being investigated. A numeri- 
cal study is posed with no such challenges and offers a 
greater degree of control in the simulation parameters. 
An arbitrary Lagrangian Eulerian approach with an 
adaptive remeshing scheme was used in [4] to simulate 
two-dimensional quasi-steady flow around liquid drop- 
lets. This method was successful in handling large de- 
formation of the interface. Three-dimensional computa- 
tions have typically found it difficult to reliably handle 
large interfacial deformations without running into sig- 
nificant numerical errors. 

In this paper, we use a moving mesh based ap- 
proach with an adaptive remeshing scheme [5]. A 3D 
unstructured finite volume scheme is used to perform 
the simulation. The interface is tracked in a Lagrangian 
manner in a deforming tetrahedral mesh. This approach 
gives us the advantage of exactly tracking the interface 
within the limits of discretization error. An adaptive 
mesh refinement scheme is implemented in order to 
maintain the quality of mesh. The adaptation involves 
mesh flipping, edge collapse and edge bisection algo- 
rithms. An optimization based smoothing algorithm is 
used in conjunction with the mesh adaptation schemes. 

Here, we present a preliminary report on the effect 
of droplet distortion on its drag coefficient. Our ap- 
proach helps us to closely correlate the deformation of a 
liquid droplet to the drag coefficient. The transient na- 
ture of the simulation will help us better understand the 
effect of flow history on the drag coefficient. At this 
juncture, we are just beginning to investigate the pa- 
rameter space and the early results provide new insight 
into the drag characteristics which have not yet been 
unearthed by experimental data. 


Governing equations and Discretization 

The governing equations for the moving and de- 
forming control volume are described by the following 
set of equations. 
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Eqn. (1) represents the conservation of volume for a 
distorting control volume. 
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Egn. (2) denotes the conservation of mass and Eqn. 
(3) gives the conservation of momentum. 
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In the above equations CV stands for the control 
volume, CS the control surface, u is the fluid velocity, v 
the velocity of the control surface, n is the unit vector 
of the face normal and f denotes the body force per unit 
mass. 

The two phases are assumed to be incompressible 
and immiscible. The continuity and jump conditions 
across the interface are given by the following equa- 
tions. 
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where o is assumed to be the constant surface ten- 
sion, R; and R> denote the radii of curvature and T is the 
tangential vector. Eqn. (4) states that the normal veloci- 
ties of both the fluids at the interface are equal and 
same as the normal velocity of the moving mesh. The 
continuity in the tangential velocity across the interface 
is given by Eqn. (5). The pressure jump balanced by the 
surface tension is denoted by Eqn. (6). The continuity 
of the viscous shear stresses across the interface is de- 
scribed by Eqn. (7). 

The governing equations (1),(2) and (3) are discre- 
tized as given in [6]. 
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where V, is the cell volume, Uy denotes the face 
flux, U oe the equivalent mesh flux of the face, A; is 


the face area, r? denotes the position of cell center 


and g is gravitational force per unit mass. Subscript n is 
the time step. 


Simulation Setup 

In the numerical investigations both the fluids are 
assumed to be incompressible, to have constant densi- 
ties py and pı, and constant viscosities, ug and for the 
gas and liquid phases respectively. The surface tension 
at the interface of the two fluids o is assumed to be a 
constant. 

The deformation factor D is defined as the ratio of 
the shortest axis to that of the longest axis of the de- 
formed droplet. The Reynolds number is based on the 
relative velocity and the gas phase properties and is 
given by, 
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The general drag force on the droplet is given by 
the Basset-Boussinesq-Oseen (BBO) equation, which is 
too complex for use in applied CFD simulations with 
large numbers of parcels. The drag force includes the 
aerodynamic steady-state drag, added mass forces, 
pressure gradient forces, and a history integral. How- 
ever, in typical CFD analyses where the droplet is con- 
sidered to be a point mass, numerous terms of this 
equation may be neglected [3]. Typically, only the 
steady state drag is used in order to define a coefficient 
of drag. 
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where 7= (;/ p is the density ratio between liquid 

and gas phases. This definition presumes that drag is 
equal to the steady state value and that the droplet may 
be represented as a sphere. The following detailed 
simulations will report the acceleration non- 
dimensionalized as in Eqn. (12). These values of Cp 
are what a typical CFD code would have to use in order 
to predict droplet acceleration correctly. These results 
are not intended to validate the correctness of Eqn. (12) 
for general use. 
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Case 1. This simulation is setup to mimic a shock 
tube experiment. A stationary liquid droplet is suddenly 
put under the influence of a gaseous flow and is impul- 
sively accelerated. The density and viscosity ratios 
maintained between the liquid phase and the gas phase 
are 50. 

The computational domain is box with a spherical 
droplet of radius r, located at the center of the box. The 
length of the box is 48r,, the width and height are 16r.. 
Boundary effects are neglected since the box is large 
compared to spherical droplet. An inlet boundary condi- 
tion is enforced on the right side of the domain and an 
exit boundary condition is set on the left side of the 
domain. 

Initially the velocity in the liquid phase is set to zero 
and a constant velocity is given to the gas phase in or- 
der to achieve a uniform velocity field. Since the veloc- 
ity is not continuous across the interface large velocity 
errors will occur as shown in Figure 1. To alleviate this 
problem the simulation is run with the droplet being 
treated as a fixed solid until the velocity jump is 
smoothed out by the viscosity as shown in Figure 2. 




















1.79 5.37 8.95 12.53 16.10 19.68 23.26 26.84 


Figure 1: Initial velocity field with stream function dis- 
continuities across the interface 


Figure 2: Velocity field smoothed out by viscosity. 


Case 2. Here a wind tunnel experiment [1] is simulated 
in which a droplet is injected in a converging wind tun- 
nel and then subjected to an accelerating flow. It is 


computationally very expensive and inefficient to per- 
form a complete direct numerical simulation of a liquid 
droplet in a wind tunnel. In this setup as well the com- 
putational domain is a box with a liquid droplet cen- 
tered in it. 

The computational domain is translated along with 
the liquid droplet in such a manner that the droplet is 
always centered. An inlet boundary condition is en- 
forced on the left side on the domain. The initial condi- 
tion imposed on the domain is similar to that of case 1. 
The same procedure to smooth velocities at the inter- 
face is followed. The velocity inlet boundary condition 
is continuously varied with time according the rules of 
flow in converging nozzles. An exit condition is pre- 
scribed to the right of the domain and the remaining 
walls are modeled as slip walls. In this case the density 
and viscosity ratios maintained are 100 between the 
liquid phase and the gas phase. 


Results 

As the droplet is impulsively accelerated by the 
gaseous flow the initial Reynolds number and Weber 
number are 40. During the acceleration of the droplet, it 
continuously deforms, the shape of the liquid droplet at 
various times is given in Figure 7. The deformation rate 
of the droplet is initially high and reduces as the relative 
velocity between the gas and the droplet decreases. The 
variation of the deformation factor with time is given in 
Figure 3. 

The Reynolds number based on the gas velocity 
goes down as the liquid phase catches up with gas 
phase. This results in the reduction on the drag coeffi- 
cient over time. Figure 4 and Figure 5 confirm this 
trend. These simulation results prove that the droplet 
deformation significantly affects the drag characteris- 
tics. 
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Figure 3: Variation of deformation factor with time. 
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Figure 4: Reynolds number vs. time 


30 





10 





Oo 25 30 
Re, 


Figure 5: Coefficient of Drag based on initial drop ra- 
dius vs. Reynolds number 
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Figure 6: Coefficient of Drag based on frontal projec- 
tion area vs. Reynolds number 


The Cp is computed on the basis of the initial drop 
radius where as C is the drag coefficient which is 


based on the frontal projection area. The volume of the 
liquid droplet will remain constant; therefore the frontal 
projection area can be calculated using the deformation 
factor data by assuming the droplet to be a prolate sphe- 
roid. It is evident that the drag predicted is much higher 
than the standard drag measured by experiments but 
correlates well with experimental data from [2]. One of 
the reasons for the higher drag coefficient is the larger 


deformation predicted by the simulation compared to 
the experiment. The numerical investigations validate 
the belief that the unsteady drag in accelerating flows is 
larger than steady drag. 

For the second case where the droplet is injected in 
a wind tunnel, the initial trends show a similar pattern 
where the predicted drag coefficients are higher than 
the ones suggested by currently accepted drag models. 
This needs to be firmly established by running the 
simulation over a larger time frame. 


Conclusions 

Numerical investigations were performed to study 
the drag characteristics of deforming droplets in accel- 
erating flows. The increased frontal area due to the de- 
formation causes an increase in the drag coefficient 
which is significantly higher than the values proposed 
by conventional drag models. These preliminary simu- 
lation results confirm this upward trend in the predic- 
tion of droplet drag. Further investigations are needed 
to study a larger range of flow conditions and the re- 
lated history effects. 
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Figure 6: Droplet shapes at times: 0.0, 2.5, 5.0, 7.0, 10.0, 12.0 and 15.0 


